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Abstract 



A microscopic theory of linear response based on the Vlasov equation is ex- 
tended to systems having spheroidal equilibrium shape. The solution of the 
linearized Vlasov equation, which gives a semiclassical version of the random 
phase approximation, is studied for electrons moving in a deformed equilib- 
rium mean field. The deformed field has been approximated by a cavity of 
spheroidal shape, both prolate and oblate. Contrary to spherical systems, 
there is now a coupling among excitations of different multipolarity induced 
by the interaction among constituents. Explicit calculations are performed for 
the dipole response of deformed clusters of different size. In all cases studied 
here the photoabsorption strength for prolate clusters always displays a typical 
double-peaked structure. For oblate clusters we find that the high-frequency 
component of the plasmon doublet can get fragmented in the medium size 
region (N ~ 250). This fragmentation is related to the presence of two kinds 
of three-dimensional electron orbits in oblate cavities. The possible scaling of 
our semiclassical equations with the valence electron number and density is 
investigated. 
PACS: 36.40. Gk, 36.40.Mr 



I. INTRODUCTION 

Complex many-fermion systems like atomic clusters and nuclei exhibit both spherical 
and deformed equilibrium shapes. Deformation of the ground state gives rise to observable 
effects in the excitation spectrum of these systems, with the splitting of the giant dipole 
resonance in deformed nuclei most likely being the best well known feature 0. Based on the 
close analogy that exists between the nuclear giant dipole resonance and the cluster surface 
plasmon, a similar splitting is expected and has indeed been observed in atomic clusters (see 
Refs. @ and || for reviews of experimental and theoretical work). A spheroidal deformation 
(prolate and oblate) is often sufficient to explain data in both nuclei and atomic clusters 
although recent jellium model calculations j|] also suggest more complicated shapes for the 
latter. 
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Clemenger || has used a deformed oscillator model, which is inspired by the Nilsson 
model of nuclear physics ||, to describe deformed clusters. However Strutinsky et al. have 
pointed out that the oscillator potential is somewhat special and argued that a spheroidal 
cavity would give a more realistic description of the equilibrium mean field (in nuclei, but 
the same is true for large atomic clusters). Actually the spheroidal cavity is still a rather 
special choice since its sharp surface is an undesirable feature. A mean field with a diffuse 
surface would be more realistic. Studies by Arvieu et al. M of the classical motion of a 
particle in a deformed Saxon- Woods like potential indicate that part of the classical phase 
space becomes chaotic. In a spheroidal cavity instead the single-particle motion is integrable 
and this fact makes calculations practicable. 

In this paper we study surface plasmons in deformed clusters by assuming a spheroidal 
cavity model for the equilibrium mean field. Thus we limit our analysis to the simplest 
deviation from spherical shape and consider only spheroidal geometry. However our model 
is not limited to small deformations and, in principle, we can study systems ranging from 
spherical to almost cigar or disk-like shapes. Our aim is to extend the semiclassical theory of 
linear response based on the Vlasov equation of Refs. PJT0[] to deformed systems. This theory 
can be viewed as a semiclassical version of the random phase approximaton (RPA). Although 
a fully quantum treatment of the problem is certainly more rigorous, the Vlasov equation 
has the advantage that the numerical effort required is greatly reduced. Indeed, a fully 
quantum RPA calculation for finite systems must rely heavily on numerical computation. 
This fact has limited explicit calculations mostly to either infinite homogeneous or spherical 
systems, where the symmetries of the mean field Hamiltonian allow for simplifications of 
the numerical problem. Pioneering work on the quantum response of spheroidal clusters has 
been made by Ekardt and Penzar |l l] (for prolate clusters only) . 

The present work is organized as follows. In Sect. II we extend the formalism of 
Refs. HlKl to spheroidal systems. In order to present clearly the main points of the theory 
we have relegated the discussion of many detailed expressions to the Appendix. In Sect. 
Ill the model is applied to study the evolution of the peak profile of surface plasmons with 
deformation. Both the "single-particle" and collective responses are studied for prolate and 
oblate geometries. Finally, Sect. IV contains a brief summary and the conclusions. 

II. FORMALISM 

In spite of its well known difficulties in reproducing the observed width of collective 
resonances, the RPA is still the basic microscopic theory of small-amplitude vibrations for 
many fermions systems. The RPA theory can be derived in the Green function formalism 
fl2]l . The RPA Green function G obeys an integral equation that can be written formally as 

G{u) = G°(u) + G°(uj)VG(u) , (2.1) 

where the unperturbed particle-hole Green function G° describes the single-particle motion 
of the constituents in the equilibrium mean field, and V is the effective two-body residual 
interaction. The response of the many-body system to a weak periodic external field of 
frequency u is proportional to the imaginary part of G(u). In practice, the calculation of G 
from Eq. ( |2.1| ) means facing two main problems: first to evaluate G° for the system under 
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study and then solving the corresponding integral equation to determine G. The latter 
becomes quite a challenge when the system is non-spherical since G will be determined by 
a system of coupled integral equations, as we shall discuss later on. 

The problem of calculating G° for a given mean field is simpler in the semiclassical 
theory of linear response developed in Refs. P,|T0fl. This theory has been applied to the 
study of giant resonances in spherical nuclei and of surface plasmon excitations in 
spherical microclusters [|14|| . The excitation spectra given by this semiclassical theory are 



very similar to those yielded by the fully quantum RPA theory even when the many-body 
system under investigation is not particularly large. Moreover for large deformations it is 
expected (Ref. 0, p. 591) that shell effects should be smaller than for spherical systems, 
thus favouring a semiclassical approach. 

The semiclassical theory of linear response based on the Vlasov equation of Refs. PJTO 



leads to an expression for the propagator with the same structure as the quantum RPA given 
by Eq. (5TT ) . However the RPA equation is actually more complicated when exchange (Fock) 



terms are properly taken into account |L2 |. Nevertheless, these terms are often treated in a 



local approximation leading to the same Hartee-like structure as the Vlasov equation. In the 
classical limit, G and G° will be denoted by D and D° respectively in order to distinguish 
them from the corresponding quantum quantities. 

Although the theory of Refs. P,|10f is generally valid for all integrable mean field Hamil- 
tonians, practical applications have been limited to spherical systems. In order to extend 
its range of application to deformed systems we consider the classical limit of Eq. fl2.1|) in 
momentum space, 

D(q', q, u) = D°(q!, q, u) + ^3 / dk D°(q', k, u) V(k) D(k, q, u) . (2.2) 



Following |9,l0], and using units % — c — 1, the propagator D° can be written in terms of 



action-angle variables {I, $} as 

D°(q[, q, u) = (2vr) 3 £ / dl F'(h (T)) 7TV~~7^~jT~- S QM > ^ Gn(q, ^ > ( 2 ' 3 ) 



n • u(I) — (u + is) 

with e — > 0, and the Fourier coefficients 



Q„(q,I) = ^ jd$e-™*e^, (2.4) 



taking the place of the quantum matrix elements (see also [IS]). In Eq. ( |2.3| ), n is a three- 
dimensional vector with integer components and the sum extends to all possible values of n. 
The components of the vector u) are the fundamental frequencies of the multiply-periodic 
particle motion in the equilibrium mean field, 

w(I) = Vjfco(I) > (2-5) 

with ho = E the equilibrium Hamiltonian. The function F'(E) in Eq. ( |2.3| ) is the derivative 
of 

F{E) = l2^ e{E r- E) < < 2 ' 6) 
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which describes the equilibrium distribution of electrons at zero temperature (Ep is the 
Fermi energy). Thus, 

f ' (£ » = -(W' (Ef - £) ' (27) 



expression that reduces by one the number of integrals in Eq. ( ]2.3|) . We also note that the 
propagator in Eq. ( |2.3| ) is perfectly well behaved when n ■ u — > due to the presence of 
the n • oj factor in the numerator. Then D° is not affected by the problem of small divisors 
(Ref. 0, p.523). 

For a realistic mean field the frequencies u do depend upon the value of the action 
variables I. This dependence reflects the non-linearity of the mean field. A cavity is an 
example where this non-linearity effect is present. Instead, for the oscillator potential model 
where the equilibrium Hamiltonian is given by 

( h ) osc = u I , (2.8) 

the frequencies & do not depend on I. This difference is a basic distinction between the 
oscillator and more realistic models and it becomes the main reason for our choice to study 
the cavity model. 

The most general partial-wave expansion of the propagator D°(q',q, uj) is 

D°(q',q,a;) = E E D LM>L , M ,(q', q,u) Y VM ,(c()YZ M (q) . (2.9) 

LM L'M' 

The symmetry properties of the system allow for a simplification of this expression. Indeed, 
the axial symmetry of spheroidal (and also spherical) systems implies 

D°lm,vm> = S m , m ,D° L 'LmW, Q, w) , (2-10) 

and therefore 

oo L oo 

£>°(q',q.") = £ £ £ Dl ILM ( q ',q,cu)Y L , M (^)Yl M (q). (2.11) 

L=0 M=-L L'=\M\ 



A similar expansion for the propagator D(q',q,u) is obtained from Eq. ( |2.2fl . Thus, for 
spheroidal systems, the three-dimensional integral equation (Eq. Q2.2| )) reduces to the follow- 
ing system of coupled one-dimensional integral equations for the coefficients D L i LM (q f , q, u>), 



1 00 POO 

— E / dkk 2 D° mM (q',k^)V(k)D eLM (k,q,u). 
j7r J e=\M\ Jo 

(2.12) 

This equation can be further simplified for spherical symmetry. In this case, 

D L , LM (q',q,u) = S L , L/ D° L (q',q,uj) (any M) , (2.13) 
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and Eq. ( |2.12|) reduces to a single uncoupled one-dimensional integral equation for each 
partial-wave component, 

D L (q\q,u) = D° L (q',q,uj) + — ^ dk k 2 D° L (q' , k, u) V(k)D L (k, q, u) . (2.14) 

The solution of this equation for surface plasmons in spherical clusters has been studied in 
Ref. [PJ], where the coefficients D L (q' ,q,u) have been explicitly derived (see also Ref. [CT). 

In order to solve the system of coupled integral equations expressed by Eq. ( |2.12| ), we 
must derive an explicit expression for the coefficients D L , LM (q',q,u). This can be done in 
a way similar to that followed in Ref. |J. We refer here to the Appendix for details on how 
to extend that approach to both prolate and oblate geometries. From Eqs. Q2.3Q and (p. 11 ) 
we obtain 



D° LILM (q',q,uj) = - 



E(2-) 3 / 



d\ 7 I de 



d(I„ 



d(E,e) 



n u u u + n v uj v + Mu v - (u + ie) 



x Q i n u ^*Mi E ^ e ^ q') Q ( n u Zl,M( E > 6 ' ^ ^) 



(LM) 



(2.15) 



where we have used the constants of motion {E, e, X z } instead of the action variables I 
to evaluate the corresponding integrals. All vector quantities in Eqs. (|2.3| ) and ( 2.15|) are 
expressed through their (u, v, <p) components as discussed in the Appendix. Furthermore, 
the integrand is evaluated at the Fermi energy Ep as a consequence of Eq. (|2.7|) . All other 
integrals are extended to the classically available phase space, the integration limits are 
specified in the Appendix. 

by 



The Fourier coefficients Q^n v m( e , e , A 2 ; q) are related to those of Eq. 

Q n (q,I)=£Q( iM V; M (q), 



(2.16) 



LAI 



and 



Qt M nln 9 = \tM,n v fdufdv 



d(u, v) 



xi L Y LM {6{u,v),0)j L {qr{u,v)) . 



[2.17) 



The symbol 



variable, the angle variable $^ takes the form $^ 



means integration over a whole period of classical motion in the respective 

&<p(u, v) + ip and ji is the spherical 
Bessel function of order L. 

The comparison of Eqs.( p.l2|) and ( 2.14T) explicitly shows that, contrary to the spherical 
case, there is a coupling between excitations of different multipolarity in deformed systems. 
It can be shown (see Appendix) that D° ULM = unless (— ) L = (— ) L , so that only multipoles 
with the same parity are mixed in Eq. ( |2.12| ). Since this is a consequence of the reflection 
symmetry of the spheroid (z — > —z), it holds both for prolate and oblate shapes. Such a 
simplification would not necessarily occur for "pear-shaped" systems. 

The physical observables we consider in this work are related to the " forward" propagator 
D(q, q, u>) regardless of the need for the full off-diagonal propagator D(q',q, uj) to solve 
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Eq. ( |2.2j ). Moreover, for randomly oriented clusters, we must average over the solid angle 
IJ. Thus we define the angle-averaged propagator 

< D{q,u) >= ^/ dqD(q,q,u;) . (2.18) 
The expansion of D(q', q, u) analogous to Eq. (|2.11|) implies 

< D(q, oj)>=^Y, D L LM(q, q, u) . (2.19) 

A quantity often used in cluster physics is the dynamic polarizability a(uS). It is related 
to the longwavelength limit of the L = 1 component of < D(q,u) >. For deformed clusters 
the polarizability is a tensor since the induced dipole moment depends on the orientation of 
the cluster relative to the external field. We introduce the M-component polarizability as 

a M {ui) = - l e hm|- [m? j . (2.20) 

In the spherical limit olm is M-independent. For spheroidal systems ck_m = olm (this is a 
property of the coefficients -Dz/la/)- 

We introduce also the "partial" photoabsorption cross sections 

47T 

omM = — wlm [a M (uj) ] , (2.21) 



and the total photoabsorption cross section is [JIT 

i 



a H = l E ^mH = | [t7o(w) + 2(7i(w)] . (2.22) 

M=-l 



III. RESULTS 

In this section we apply the formalism developed in Sect. II and in the Appendix to the 
study of prolate and oblate clusters with varying degree of deformation. Also we explore 
the possibility of scaling with the number of valence electrons and cluster density. The 
emphasis is on understanding how the collective response is altered when the clusters under 
investigation do not have a spherical shape. We shall not attempt to compare our results 
with experimental data at this level of development. For this reason we do not consider 
such refinements as the electron "spill out" which is well known to give a red shift of the 
plasmon peak in spherical clusters. We expect to find a similar effect in the deformed case. 
No particular effort is made to reproduce the observed width of the plasmon resonances 
which is considerably larger than that obtained in the present model. 

The picture we have in mind here corresponds to taking a spherical cluster and then 
deforming it to either a prolate or oblate shape while keeping its volume and density constant. 
We consider initially a relatively large sodium cluster of spherical shape, with iV = 254 
valence electrons, described approximately by a square-well mean-field potential of radius 
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R = r s and r s is the Wigner-Seitz parameter in units of the Bohr radius. We will take 
r s = 4.0. In this case the surface plasmon resonance consists of a single peak situated near 
the Mie frequency, 

^Mie = ^= , (3.1) 

with uj p the bulk plasmon frequency. Then, this sodium cluster is deformed to a spheroidal 
(prolate or oblate) shape characterized by a deformation parameter 77, 

V=j£, (3-2) 

where f?> (R<) are the larger (smaller) diameters of the spheroids. We study the changes 
that are expected in the distribution of the dipole strength and report results for 1 < 77 < 2. 
The spherical limit is obtained by taking 77 = 1.001. 

To describe the spheroidal systems we introduce spheroidal coordinates {u, v, f} as de- 
scribed in the Appendix for the prolate and oblate geometries. Assuming that electrons 
move in a static cavity of spheroidal shape, the equilibrium Hamiltonian /i will be chosen 
to be 

h (u, v, (p, p u , p v , p v ) = K(u, v, <p, p u , p v ,p (p ) + V(v) , (3.3) 

where {p u ,Pv,Ptp} are the corresponding conjugate momenta, K is the kinetic energy and V 
the potential energy for particles in a cavity. All the calculations have been performed within 
this framework. The central objective has been to calculate the collective Dum function 
from Eq. ( |2.12| ). The number of coupled integral equations required is 3 for the cases of 



larger deformation (77 = 2) considered in this work (that is i = 1, 3, 5 in Eq. ( |2.12j )). 

One technical point deserving special attention is the handling of the singularities present 
in the D° propagator in Eq. ( |2.15 ). We have taken a finite value e = 0.0020^6 when per- 



forming the actual calculations. Our choice modifies slightly (~ 5%) the width of either the 
single-particle or the collective plasmon peak. Furthermore, we have summed up all allowed 
frequency modes for — 2 < {n u ,n v } < 2 when calculating D°. The dynamic polarizability 
and the total photoabsorption cross section of prolate clusters are mainly dominated by the 
n u = 1, n v — 0, ±1 modes for M = and n u = 0, n v = 0, ±1 for M = 1. However, the 
oblate cavity presents some peculiar features requiring a finer analysis. 

From a theoretical point of view it is interesting to study also the single-particle response 
(proportional to the imaginary part of the zero-order propagator D°) since its features are 
more directly related to the shape of the equilibrium mean field. This can be of help in 
understanding how this shape affects the collective response. 

The effective two-body residual interaction V(k) in Eq.( |2.2| ) determines primarily the 
position of the collective plasmon peak and, to a lesser extent, its shape. The present 
semiclassical approach is essentially a Hartree approximation which does not include ex- 
change contributions. However, exchange and correlation terms can be taken into account 
in a local approximation by introducing a (static) local-field correction G(k) |T7|]. Thus the 
momentum-space interaction we use is 

V(k) = 4n^ [1-G{k)\ , (3.4) 
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with 



G(k) = A 



(3.5) 



where A = 0.9959 and B = 0.2612 for sodium ( ||17|| ) p. 446) and pf is the Fermi momentum. 
In our calculations we have verified that it is possible to take an upper integration limit 
kmax — ^Vf when solving Eq.( [2.12| ) without a noticeable change in the relevant results. 



A. Prolate cavity 

The natural frequencies {u u , u> v , u v } of the unperturbed trajectories of particles in the 
prolate cavity determine the gross behaviour of the semi classical propagator D° (Eq. (|2.15| )) 
as a function of uo. Since these frequencies depend on the integration variables e and X z it is 
interesting to evaluate the possible values that they take and their occurrence as a function 
of the cluster deformation r). In Fig. [TJ we show histograms for the probability density of 
the natural frequencies expressed in terms of the Mie frequency oJuie and for deformations 
7] = 1.0, 1.5 and 2. This probability density refers to the occurrence of the frequencies as 
functions of the (discretized) integration parameters e and A 2 in our numerical calculations. 
It was evaluated with bins of width 0.002 uiMie- The frequencies uj u and uj^ fall in a narrow 
range of the order of 0.05 uJuie-, with the former moving to the left and the latter to the right 
of the spherical limit as deformation increases. For uj v we have not plotted the symmetric 
negative components in Fig. [I]. The frequency u v shows a different behaviour and it spans, 
in principle, the [~ 0.5 uiMie, oo] range. Overall, Fig. [I] gives a clear idea of the dominant 
frequency poles contributing to the D° propagator. 

The M = 0, 1 components of the photoabsorption cross section (Eq. (|2.21|) ) per valence 



electron are shown in Fig. ^| for prolate clusters with deformation ranging from rj = 1.0 to 
2., step 0.25. We display the results in a 3-dimensional plot to better assess their relative 
behaviour. The areas shaded black show the cross section obtained from the zero-order 
propagator D°. As a reference, we can observe on the left-hand side of Fig. ^| the same 
peak corresponding to spherical geometry (77 = 1.0) in both M = and 1 components and 
centered around u> ~ 0.16 oj Mie . Then, the M = component shows a photoabsorption 
peak shifted to lower frequencies as deformation increases. This is a consequence of u u (see 
Fig. [1|) determining the dominant pole. The corresponding strength of the peak increases 
for r] = 1 to ~ 1.5 and then slowly starts to decrease. The M — 1 peak is shifted to 
higher frequencies due to the dominance of uj^ in this mode and its strength decreases with 
increasing deformation. 

The right-hand side of Fig. Q shows the partial collective cross sections (shaded gray) 
evaluated from the full propagator D in Eq. ( |2.20| ). By including the fluctuation of the 
mean field we obtain two main effects: a huge shift in the position of the M component of 
the photoabsorption peak to frequencies around the Mie frequency and a noticeable change 
both in the width and the strength of the peak, mainly for M = 1. We have numerically 
checked that the energy- weighted sum rule (area under the curves) is unchanged within 2%. 
The low-energy collective plasmon (M = component) has an intrinsic structure which is 
slightly more complex than the corresponding high-energy plasmon (M = 1 component) 
and showing some degree of fragmentation noticeable at rj ~ 1.5. The dominant collective 



S 



peaks in Fig. |] are slightly blue shifted with respect to the positions predicted by the Mie 
theory both for the M = and 1 components. This completely classical theory predicts that 
the M = and M = 1 peaks should be at the frequencies u>i corresponding to oscillations 
along the i-axis [z for M = 0, x and y for M = 1). These frequencies are given by 



UiUJp , 



with rii the appropriate depolarizing factor. For spherical symmetry n a 



giving Eq. while for a prolate spheroid |T8| 

l + e 



n 



y 

1-e 

2e 3 



log 



2e 



The eccentricity e is related to our deformation parameter 77 by 



1-4 . 



(3.6) 



1 

3' 



(3.7) 
(3.8) 

(3.9) 



A simple calculation gives u> Xty = 1.073 u>Mie, = 0.836 u Mie for 77 = 1.5 and u) x>y = 
1.113 WMie, ^2 — 0.722 o;jwj e for 77 = 2. From Fig. |2] we observe that our collective plasmon 
peaks are blue shifted by ~ 7% with respect to the position expected from the classical 
values. 



B. Oblate cavity 

We turn now to present our results for oblate sodium clusters. One of the most interesting 
aspects of the present cavity model is the existence of two kinds of three-dimensional orbits, 
W and B orbits, for oblate geometry (see Appendix). This feature should be shared, at 
least qualitatively, by more realistic deformed mean fields. In order to estimate the relative 
importance of the two kinds of orbits, we have plotted in Fig. |3| the fraction of valence 
electrons following each of the orbits as a function of the cluster deformation 77. In the 
spherical limit, 77 — > 1, the cavity only allows W-type three-dimensional orbits since the 
B orbits become oscillations along a diameter. As deformation increases the fraction of 
electrons in B orbits raises quickly reaching 50% for 77 ~ 1.7. The presence of W and B 
orbits in the oblate cavity suggests that we may expect a richer behaviour of the cluster 
response. 

The W and B orbits are characterized by their respective fundamental frequencies 
(Eq. ( |2.5| )) which in turn determine the D° propagator. We shall refer to them as the 
and {u) B } frequencies. In order to have an estimate of the possible values they take and of 
their occurrence for different deformations of the oblate cluster we present histograms for the 
probability density of the natural frequencies {u^ ,lo^ ,u^} in Fig. [| and of {lo b ,u^} 
in Fig. H]. The {u^} frequencies in Fig. [| may be compared to the natural frequencies ob- 
tained for the prolate cavity shown in Fig. |l]. Clearly the results for 77 = 1 in both cases are 
identical. A great similarity is also observed for u v as the clusters get deformed. However, 
the increasing occurrence of uj v in Fig. |1] at the lower limit of the frequency range is not 
present in (Fig. |j). The allowed values for 10™ change slightly with deformation and are 
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not sharply defined. The behaviour of uo u in Fig. [I] is rather different, with the frequency 
defined in a narrower range and its allowed values decreasing with increasing deformation of 
the prolate cavity. The uu^Y frequency in Fig. ^ shows a similar structure to its counterpart 
in the prolate case but its magnitude decreases with deformation rather than increasing 
(Fig- 1). 

The {uj b } frequencies are characterized by a different probability density , as shown in 
Fig. In the spherical limit (r] — 1) the B orbits have perfectly defined frequencies and 
u B = ± oj b = | u B . As deformation increases ui B becomes much less defined and, on average, 
takes smaller values than in the spherical case. The frequency uj b instead is still fairly well 
defined and is also decreasing with increasing deformation. Furthermore, uj b shows great 
similarity with uo^ in Fig. f|. The frequency uj b shows a definitely different behaviour than 

. It spans a finite range of values depending strongly on the cluster deformation. From 
the results shown in Figs. |4| and f| we may expect a more complex behaviour of the D° 
propagator than that observed in the prolate cavity. 

The M — 0, 1 components of the photoabsorption cross section (Eq. ( [2.21| )) per va- 



lence electron are shown in Fig. |6| for oblate clusters with deformation ranging from 
r] = 1.0 (0.25) 2.. As in Fig. g we have plotted the cross sections cr^J (shaded black) ob- 
tained from the propagator D° calculated in the static mean field and om (shaded gray) 
calculated from the correlated propagator D. For M = we observe presenting a 
novel structure. Indeed, the dipole peak is splitted into several peaks as deformation of the 
oblate cluster increases. A simple estimate based on the probability densities for u w and 
u B indicates that both W (dominant) and B orbits contribute to the left most peak in the 
mode {n u — l,n v — 0}. At higher frequencies the other two peaks come from B orbits 
in modes {n u — — l,n v — 1 } and {n u — l,n v — 1} respectively. Once fluctuations of the 
mean field are taken into accout, we obtain a rather complex structure for cr (upper gray 
peaks). As deformation of the cluster sets in, the M = component of the collective dipole 
peak corresponding to the spherical cluster (77 = 1.0) gets fragmented as a consequence of 
the more complex structure of D° in the resonance region. The amount of fragmentation 
depends on the deformation parameter. The high frequency peak is more fragmented for 
1 < T] < 1.5, reaching a rather simple structure with a dominant peak for r] = 2. The M = 1 
component of the photoabsorption cross section shown in Fig. does not present major new 
features. The a\ shows a dominant peak for each deformation {n u = 0, n v = 0) with W 
and B orbits contributing in proportion to the number of valence electrons moving in each 
of them (see Fig. |3|). The collective effects shift the o\ dipole peak to frequencies below the 
one corresponding to the spherical cavity as deformation increases. 

The position of the dominant peaks can be compared with the predictions of the classical 
theory. In the oblate case, the o>j frequencies are still given by Eqs. (|3.6D (P7T|) and |]TS 



1 + e 2 r . . . 

n z = — - — e — arctane , (3.10) 

e 6 

with the eccentricity e related to our oblate deformation parameter by 



7/ 2 -l. (3.11) 

A simple calculation gives ou X) y = 0.912^^6 , w z = 1.157 ujuie for V = 1-5 and ou X) y = 
0.842 u>Mie j <^z = 1-258 c^Mie f° r V — 2. The M — 1 collective plasmons (u x ,u y ) are 
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again blue shifted and in roughly the same amount we observed in the prolate case (~ 7%). 
However for M = the collective peaks are now blue shifted by ~ 10 — 15% with respect to 
the position expected from the classical values. 



C. Comparison between prolate and oblate clusters 



In this section we compare the most prominent features concerning the dipole response 
of spheroidal clusters to an external field. We have chosen to show the photoabsorption 
cross section per valence electron (N = 254) calculated for a sodium cluster which has 
been deformed to both prolate and oblate shapes. We display calculations corresponding 
to deformations rj = 1.0(0.25)2 in both geometries. In Fig. |7] we plot the corresponding 
cross sections calculated both in the static mean field (single-particle approximation, shaded 
black) and including collective effects (shaded gray). Arrows are plotted as reference. They 
place the position of the plasmon peak for a spherical cluster in the present model. For 
prolate geometry the photoabsorption strenght always displays the characteristic splitting 
into two pronounced peaks, corresponding to oscillations along two perpendicular symmetry 
axes (M = and M = 1). The relative strength of the two collective peaks varies with 
an increasing dominance of the M = 1 component for larger deformations. For oblate 
spheroidal geometry instead, an interesting phenomenon occurs: one peak dominates the 
cross section in the rj = l-=-2 range. The high-energy collective plasmon (M = component) 
gets fragmented and distributed over an interval of frequencies of the order of ~ 0.3 uouie- 
This fragmentation tends to dissapear at the largest deformation studied here [rj = 2). 
However, as shown in Fig. some residual fragmentation remains at frequencies similar to 
the characteristic frequency of the low-energy plasmon (M = 1 component). The effect has 
been traced back to the appearance of extra peaks in the "single-particle" M = strength, 
that is related to the existence of two kinds of three-dimensional orbits (W-type and B-type) 
in an oblate cavity. It is reasonable to expect that a similar effect should occur also for more 
realistic mean fields. 

Overall our results give a clear indication on the gross features to be expected for either 
prolate or oblate medium size sodium clusters. The photoabsorption cross section is splitted 
into two peaks for the prolate geometry, but is mainly dominated by a single peak in the 
oblate case for the range of deformations studied here, the smaller peak being further splitted 
or fragmented. 



In this section we investigate the dependence of our results from the valence electron 
number N and from the density parameter r s . By using the explicit expressions of the free 
propagator D° given in Sect. II and in the Appendix it is possible to define for each geom- 
etry (prolate, oblate, spherical), at given cluster deformation rj, a corresponding universal 
propagator <i°(x / ,x, s) such that 



D. Scaling properties 




(3.12) 
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where uj p = -^g- and R is the cluster radius for spherical geometry. In Eq. (|3.12j ) all the 
dependence from the number of valence electrons and from the electron density has been 
extracted from the free propagator D° by relating it to a universal function d° of the 
dimensionless parameters x = i?q and s — — . The question we address now is whether a 
similar factorization occurs also for the collective propagator D. In general this factorization 
is not possible for an arbitrary interaction V. Indeed, by defining a function d(x.', x, s) from 
the collective propagator D through Eq. ( |3.12j ) then the integral equation ( |2.2|) leads to the 
following integral equation for d, 

d(x', x, s) = d°(x', x, s) + — J dy d°(x', y, s) V (|) d(y, x, a) , (3.13) 

which explicitly contains the density parameter r s . An additional N and r s dependence is 
implicit through the R parameter in the force V. Consequently it is not possible to determine 
an obvious scaling law allowing for the calculation of a universal collective function d(x', x, s). 

In our case V corresponds to a modified Coulomb-type residual interaction (Eq. ( |3.4| )). 
Then it is reasonable to search for additional properties of Eq. (|3.13| ) in the presence of a 
pure Coulomb force. We obtain 



d(x',x,s) =d°(x , ,x,s) + — ^r s N~3 I dyd°(x',y,s)^^d(y,yi,s) . (3.14) 



2 



(27r )3— J ^">-> y 

In this limiting case, the solution to Eq. ( |3.14| ) depends on the product r s iV§ only. De- 
viations from this result are due to more realistic effective interactions and they may be 
expected to be important in the smaller clusters. 

In the limit of large N our microscopic calculations tend to reproduce the results of the 
macroscopic theory, with the plasmon frequency depending only on the electron density, 
but not on the total electron number N. The iV- dependence of Eq.( [3.14j ) is not in conflict 
with this expectation since the dimensionless parameter s is related to the frequency by an 
iV-dependent relation s oc N$ r 2 s LU. 

Our discussion indicates that no exact scaling is to be expected for the photoabsorption 
cross sections in the present model. The actual amount of scaling violation in D does 
depend on the size, composition and geometry of the particular cluster considered. In order 
to illustrate the iV-dependence of our results we have calculated the collective response of 
sodium clusters with 25 atoms, thus changing iV by an order of magnitude with respect to 
the results reported earlier. In Fig. [8] we display the collective response of iV = 25 sodium 
clusters (shaded black) together with that for iV = 254 (shaded gray). For spherical and 
prolate geometry the N = 25 plasmons are blue shifted by ~ 10% with respect to the 
iV = 254 results, with the overall peak profile showing a qualitatively similar behaviour in 
the two cases. The light cluster plasmons are sharper due to a smaller Landau damping. 
For oblate geometry we notice that the fragmentation of the high-frequency plasmon found 
for TV = 254 does not appear for iV = 25. The fragmentation phenomenon is iV-dependent 
and it disappears for different cluster mass. These results agree with our general statement 
about the lack of scaling with iV of the collective propagator D. 

The study of clusters with a larger number of valence electrons (N > 254) will not 
provide much new information since the theory starts approaching the classical limits. We 
have confirmed this by doing calculations for sodium clusters with iV = 2500. We observe a 
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red shift with respect to the N = 254 case of the order of 3% and no fragmentation of the 
M = peak. 

We have investigated also the electron-density dependence of our results. Calculations 
of the photoabsortion cross section for prolate and oblate sodium and potassium clusters 
(r s = 5) at iV = 254 show very similar peak profiles (including fragmentation) when the 
frequency axis is rescaled with the appropriate value of lOuie- Although the density of valence 
electrons was decreased by about a factor 2, these small variations can be understood on 
the basis of Eq. (|3.14 ) where the relevant parameter (r s iVs) was changed by 25%. 



IV. SUMMARY AND CONCLUSIONS 



In this paper a semiclassical theory of linear response based on the Vlasov equation PJlQ 



has been extended to spheroidal systems and then applied to study the peak profile of surface 
plasmon resonances in medium size deformed atomic clusters. Assuming a spheroidal cavity 
model to describe both prolate and oblate clusters we have been able to calculate the gross 
features of the cluster response to an external field of frequency wasa function of the cluster 
deformation. 

Two main general results emerge from our calculations. One of them refers to the 
splitting of the collective dipole peak with increasing deformation and to the position of 
these peaks. On this we have commented at lenght in Sec. III. The other main result is 
related to the width of the dipole peaks. In the present model the single-particle dipole 
resonance does display a width which is due to the non linearity of the assumed equilibrium 
mean field. This single-particle width generates a width in the collective plasmon resonances 
through a mechanism that is analogous to the Landau damping in homogeneous systems. 
Our calculated width is not sufficient to reproduce the observed plasmon width. However 
our width is underestimated since we have included only the first few frequency modes when 
evaluating D° in Eq. ( p.3|) . The neglected terms would increase the imaginary part of D° 
in the region of u ~ ojMie and thus increasing the Landau damping. Estimates based on 
numerical calculations including more modes set this effect at about 20%, which is far from 
sufficient to explain the observed values. Thus more sophisticated effects, like the possible 



coupling to surface vibrations |jT9| , should perhaps be taken into account. 

The fragmentation of the high-energy plasmon peak for medium size oblate clusters is 
a definite prediction of the present theory and reflects the existence of nontrivial dynamics. 
The relative importance of this effect depends on several elements such as the shape and 
size of the cluster and the strength of the effective interaction. 

The normalized cluster dipole response shows no simple scaling properties with N at 
given density (r s = 4). The position of the plasmon peaks is weakly dependent on N, 
changing by ~ 10% in the 25 < iV < 2500 mass region while fragmentation of the high 
energy component dissapears for light (N ~ 25) and heavy (N ~ 2500) oblate clusters. On 
the other hand, no major changes in the photoabsorption profile are observed when density 
decreases by a factor ~ 2 (r s = 5) for a given number N of valence electrons. 

Extensions to this work are clearly welcomed. In particular, a close comparison to ex- 
periment is required to asses in detail some of the physics missing in the present model. 
Nevertheless we have established a solid framework to classically understand the gross fea- 
tures of complicated, intrinsically quantum systems. 
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APPENDIX: 

In this Appendix we specify the details for the expressions given in Sect. II. The propa- 
gator D° is determined entirely by the single-particle motion in the equilibrium mean field. 
The classical motion of a point particle in a spheroidal cavity with perfectly reflecting walls 
has been studied in Ref. J7J for prolate cavities and in Ref. || for both prolate and oblate 
cavities (see also Refs. [pOjpI] ). The authors of Ref. [§] have pointed out that the three- 
dimensional motion in a prolate cavity is simpler than that in an oblate cavity. We give 
first a detailed description for prolate geometry. Then, through a simple transformation, 
the corresponding results for the oblate shape can be recovered. Furthermore we remark 
the non trivial differences related to particle motion in the two geometries. Our notation 
follows closely that of Ref. |7| . 



1. Prolate cavity 

To describe the prolate shape, let us introduce the prolate spheroidal coordinates {u, v, f} 
through their relation to the cartesian coordinates, 

x = £ p cos u sinh v cos ip , — | < u < f , 
V — C P cos u si n h v srn f > < i> < oo , 

z = £ p sin u cosh v , < (p < 2n , (Al) 

with 

£ p = v^3^2 _ ( A2 ) 

In the last equation c and a are the larger and smaller semiaxes respectively. The two focal 
points are at z = ± £ p . 

The equilibrium Hamiltonian ho can be expressed in terms of the spheroidal coordinates 
{u,v, ip} and their conjugate momenta {p u ,Pv,P<p} EH- Assuming that particles move in a 
static cavity of spheroidal shape, represented by a potential energy V(v), then h is given 
by 

2 2 2 

h ( p \u,v, <p, Pu ,p v , P(p ) = - P « + P " , 2 + - .% — + V(v) , (A3) 

2 m (cosh v — sin u) zm^ z p sinh v cos z u 

for the prolate configuration. The potential energy for sharp walls is 

V{v) = v < v% , 

= oo v > t»f , (A4) 
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and the v 2 parameter is determined by the shape of the cavity, 



sinh-uf = — . 



(A5) 



The Hamiltonian (Eq. ( |A3|) ) is integrable and the particle motion could, in principle, 
be described in terms of the angle variables {$ u , $„, and of the three conjugate action 
integrals 



-!- / p u du 



1 

2^ 



p v dv 



h fa*** 



(A6) 



However, in the spirit of Ref. |J, the following three other constants of the motion {E, e, X z } 
can be conveniently used instead of the action integrals. These new constants are: the 
particle energy E, the separation variable e and the z-component of the particle angular 
momentum X z (which coincides with the action variable I v and with the generalized mo- 
mentum ptp). The constant of motion e plays a role analogous to the magnitude of the 
particle angular momentum in the spherical case. 

With the help of the Vlasov equation the angle variables can be explicitly expressed in 
terms of these three constants of the motion {E, e, X z } and of the spheroidal coordinates 
{u, v, f}. The derivation, based on separation of variables in the linearized Vlasov equation, 
is lengthy but straightforward, hence we do not report all the details here W2 . 

The generalized momenta p u v are 



Pu 
Pv 



£ p V2^Ej<rf-U p (u,a£) 



i p y/2^E^/V p (v,ai)-a[. 
The two dimensionless constants of the motion erf and ovf are defined as in Ref. 



e 
~E ' 



and the effective potentials in Eq. 

U p (u, erf) = sin 2 u + 
V p (v, erf) = cosh 2 v ■ 



2 2mEQ ' 
for the u and v coordinates are given by 



(A7) 



(A8) 



cri 



COS 2 u 



af 



sinh v 



(A9) 



Particles on three-dimensional orbits move between two confocal ellipsoids with vf < 
v < v£ and uf < u < uf , with v[ 2 an d uf 2 the turning points obtained as solutions to 
p v = and p u = respectively (see Ref. for details). Introducing the quantities 



\ 



cri 



± 



\ 



1 



(AlO) 



these turning points are obtained from 
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COsh«) = y/l- {t P _f , 

cos(m^) = t+ , % = — u 2 , (All) 

and from Eq. ( |A5|) . 

It is our purpose here to specify the integration limits in Eq. Q2.15|) , and to derive explicit 



expressions for the eigenfrequencies uJ u ,v,<p as well as for the Fourier coefficients Q^$ v ,n 
appearing in the expression for D°. 

To obtain the integration range for the constant of motion e we realize that the turning 
point U2 exists only if cosii^ < 1 in Eq. ( All| ). Also the presence of an infinite poten- 



tial barrier at the surface of the cavity implies p„(i>2 ) — 0- These two conditions impose 



constraints on erf, 



°1 ) ■ = a 2 > 

' mm 



of) =cosh% 2 p -^H,. (A12) 
) max sinh 1>2 



These constraints determine the integration range for e. Thus, 



de — > de, (A13) 



with 

- = 2^f <+ = ^ C0Sh2 <-^lv- (AM) 

The integration range for the constant of motion A z is determined from the possible 
values of the particle angular momentum along the symmetry axis. Thus, 

/ d\ z — > £ + d\ z , (A15) 

where 

A± = ± ( PF a) , (A16) 

with pf the Fermi momentum associated to the Fermi energy Ep. 

The evaluation of the Fourier coefficients in Eq. Q2.17 ) can be simplified by noticing that 



n u ® u + n v $ v + n v $ v = s n (u) + s' n (v) . (A17) 
The phases s n (u) and s' n (v) are given by 

s n (u) = n • u t u (u) + n ■ Q a u (u) - n v ^/ u (u) , 

s' n (v) = n • wt v (v) + n • Qa v (v) - n^^v) . (A18) 

with u> the frequency vector defined in Eq. (|2.5|) , Q, an auxiliary frequency vector and the 
auxiliary functions { a uv , t U)V , 7 Uj „ } given by 
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a u (u) = m£ 



du' 
Pu{u') 



. . , 2 f u sin 2 ?/ 

t u {u) = - rn£ p / - 

Ju m i n Pu\ U ) 



-du', 



~f u (u) = X z 



du' 



u mm COS 2 U'p u (u') 



a v (v) = — m£ 
t v {v) =m£l 
lv{v) = X 



dv' 



Vmin Pv{v ') 

v cosh 2 t> ' , . 

— r^r ' 

rn i ri 

dv' 

sinh 2 v 'p v (v ') 



(A19) 



In these equations the upper integration limit is a variable while the lower one is the corre- 
sponding inner turning point for the geometry under study. For the prolate case u min = uf 
and v min = t?f . All the integrals appearing in Eq. ( |A19| ) can be easily expressed in terms of 
elliptic integrals and evaluated numerically 

The natural frequencies of the unperturbed trajectories are the three components of the 
frequency vector u (Eq. ( |2.5| )). They can be expressed in terms of these integrals as 



u u (E, e, X z ) = n 



UJ V (E,€, \ z ) = 7T 



a u {u^)r v (v^) - a v {v^)T u (u^)} ' 



7T 



uj v (E, e, \ z ) = — u)vlv{v%) +a; u 7„(wf) 



(A20) 



Correspondingly, the auxiliary frequency vector f2 has components 



VL U (E, e.X z 
Q V (E, e.X z 



— U), 



t u (u%) 



a u (u%) 



TV 



(A21) 



Note that from the expressions given above we can easily obtain the explicit form of 
the angle variables in terms of the spheroidal coordinates. For example the explicit form 
of v) can be derived from Eq. ( |A17[ ) by setting n u — 1, n v — 0, n v = 0. The angle 
variable takes the form $ 9 = $ip(u, v) + tp. 

The evaluation of the Fourier coefficients in Eq. (|2.17| ) is made explicit by using the 
identity 



du j> dv 



du 



dv 



d(u, v] 



- i (n u $u +n v 3>„ +n v <E> v ) 



d(u, v) 



cos[s n (u)] cos[s^(u)] , 



(A22) 



and evaluating the Jacobian, 



d(u, v) 



2 U!., U),; 



( C 2 Y M u 



PuPv 



a u (u$) a v (v%) 



(cosh 2 - ^- 2 



v — sin u 



(A23) 



17 



Then, 



f u 2 du f v 2 dv 

Q { n u f v ,n v = hl,nj p — — F(u, v) COs[s n («)] COs[s' n (v)] 



(A24) 



with 



F(u,v) 



,r 4 / \2 



cosh 2 v — sin 2 u 



x y £A/ (0(u, u), 0) j L (qr(u, v)) . 



(A25) 



The radial coordinate r and the polar angle # can be easily expressed in terms of the (u, v) 
variables by using the relations (|Al|). 

The present calculation is simpler if parity selection rules are taken into account. These 
selection rules originate from the fact that the effective potential U p (u, af) is an even func- 
tion of u and as a direct consequence of the spheroidal geometry invariance under the 
reflection z — > —z. Because of this symmetry some of the Fourier coefficients (Eq. ( |2.17| )) 
vanish. This can be proved by using the following relations, 



Sn{-u) = n u ir - s n {u) , 
cos[s n (-*u)] = {-) nu cos[s n (u)] , 
cos[6 l (— u, v )] = — cos[9(u, v )] = cos[7r — 9(u, v)] 
Y LM (e(-u,v),0) = (-) L - M Y LM (6(u,v),0) . 

All the remaining factors in (|A24j) are even functions of u. Then 



(A26) 



rutin dii r v o fin 

Qt M ln v = WJ1 + (-r +L ~ M ] / - f ^F(u,v)cos[s n (u)]cos[s'M], (A27) 

Ju p u Jvf p v 



with uq = 0. Equation ( [A 2 7] ) implies that, for example, for L = 1 and M = 1 we need to 
sum only over even values of n u , while for M = we must take only odd values of n u . It 
means also that D° ULM = unless (— ) L = (— ) L . 



2. Oblate cavity 



The transformation analogous to ( |A1| ) for oblate geometry can be obtained from those 
equations by exploiting the identities 



sinh(f ± i | ) = ± % cosh(t> ) , cosh(t> ± i |) = ± i sinh(i>) . 



(A28) 



Then the following formal replacements should be made in (|A1 ) to obtain now the relation 
between the oblate spheroidal coordinates {u,v, tp} and the cartesian coordinates, 



it 



O ' 



v-i%, 



(A29) 



where 



Co = Va 2 - c 2 



(A30) 
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is the radius of the focal circle and a is now the largest semiaxis. Applying the same 



transformation to the kinetic enegy part in Eq. (A3) we obtain the corresponding equilibrium 



Hamiltonian for an oblate cavity. The potential energy is still given by Eq. ( |A4| ) , with 
being replaced by v 2 °, 



cosh-u 2 ° = — . 



(A31) 



Note that this relation can be obtained by applying the transformation ( |A29| ) to Eq. ( A5). 

The generalized momenta p U)V are obtained by applying the transformation ( |A29|) to the 
expressions in ( |A7| ) and by making the replacements 



.o 



of, 



a. 



o 



2m EQ 



a p 2 



(A32) 



The generalized momenta for the oblate cavity are then expressed as 



Pu 



Pv 



i V2rnEJa°-U (u,a°) 



£ V2m~EjV c 



(A33) 



with the oblate effective potentials 



U {u,a 2 



V (v,a 2 °) 



2 ^ *2° 
cos u + 



COS 2 u 



cosh v + 



a. 



o 



cosh v 



(A34) 



Of course we still have = A 2 



In Ref. H it has been pointed out that the effective potentials U Q (u, a 2 °) and V Q (v, cr 



for oblate geometry may exhibit a non monotonic behaviour for some range of values of 
the parameter a 2 °. As a consequence the phase space of an oblate spheroidal cavity is 
divided into two parts and there are two kinds of three-dimensional orbits (plus a separatrix 
that, however, has zero weight in our calculations). For their description we introduce the 
equivalent to the quantities t± in Eq. ( |A10|) for the oblate case, 



■ o 



\ 




(A35) 



which are convenient for expressing the turning points. Of course these parameters can be 
obtained by making the replacement ( |A32| ) in Eq. (|A10 ). 

We must now distinguish between the two kinds of three-dimensional orbits occurring 
in the oblate cavity: 



B— orbits 

For 



a? <1 



and 



'<7 2 ° < O 



o 



< 1 



a 



o 



(A36) 
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the orbits always cross the focal circle. These are the orbits with a hyperboloidal 
caustic of Ref. || and they are analogous to the so-called bouncing ball modes. In 
this case, the equation p u = has four solutions. The accessible region of phase space 
is uf < u < and uf < u < u 2 B and with < v < v 2 °. For the lower integration 
limit in the v variable we have to take vf = since the equation p v = has no real 
solution in the interval [0, v 2 ]■ The explicit turning points for the u variable are 
obtained from 



cos U 2 = t° , 



cos-u. 



, B 



t 



o 

+ ' 



B 



— Ur, 



!/■! 



. B 



- U A 



. B 



(A37) 



The constraints expressed in Eq. ( |A36| ) imply the following integration range for X z 
and e, 



d\ 7 



de 



de 



(A38) 



with 



± (PF £ 

2 Ei 



O) ' 



E F (1 + o, 



(A39) 



W— orbits 

For 



a? > 



and 



a 



o 



> 1 + a. 



o 

2 5 



(A40) 



the orbits have an ellipsoidal caustic, and they are analogous to the so-called whispering 
gallery modes. For these orbits p\ > for < u < and pi > for v± < v < v! P, 
with the corresponding turning points given by 



cos-u^ = t 



2 
W 



U i 



. w 



Un 



, w 



cosh Uj = t 



■ o 



(A41) 



The integration limits for the constant of motion integrals ( cf. Eq. (|A38|) ) are 

= ± (p F a) , 



. w 



E F {l + o°) 



. w 



E h 



. o (To 

cosh v + 



o 



cosh V2 



o 



(A42) 



The integration limit is determined by the cavity surface. 
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The number of particles moving on each kind of orbits can be easily evaluated with an 
integration over phase space. The total number of valence electrons is 



with F(E) given by Eq. 



N 



N = J d$ dl F(E), 
Clearly, setting U — 1 
h (I)) , 



(A43) 



2 / dl6{E F 







dE I dX z I de 



d(I v ,I u ) 



d(E,e) 



(A44) 



We can define N = N B + N w with N B (N W ) the number of valence electrons in i?(W)-orbits 
respectively and 



dE 



w 



dE 











d\ z 






'e1(E) 










d\ z 









de 



de 



d{I v ,I u ) 



d(E,e) 
d(I v ,I u . 



d(E, e) 



(A45) 



The limits for the A and e integrals in these formulae are given by Eqs. ( |A39|) and (|A42|) 
where Ep is replaced by E. Since the propagator ( |2.3| ) has the same structure of an integral 
over the classical phase space as N, it is convenient to make a similar distinction between 
the contribution of the two kinds of orbits. Thus for oblate geometry we write 



d l'lm(Q^ w ) = B L r LM (q', q, u) + W ULM (q\ q, u) , 



(A46) 



with Bl'lm and Wl'lm still given by Eq. (|2.15|) but with the corresponding integration 
limits for each kind of orbit. 

The building blocks of the present calculation are the elliptic integrals given by Eq. (|A19|) 
for the prolate case. They determine the frequencies ( |A20| - |A~2~T| ) , as well as the phases ( |A18| ), 
required in the evaluation of the Fourier coefficients. The analogous expressions for oblate 
geometry are obtained by applying the transformation (|A29|) to the prolate formulae. All 



formulae given for prolate geometry can be translated in the same way for the oblate case. A 
little extra care should be taken for the u integrals in the case of B-orbits since the double- 
well structure of the effective potential U Q (u, <t 2 °) makes p u become imaginary in the interval 
[u 3 , u^]. This integration range must be excluded by very definition of classical phase space. 



Then, the lower integration limit for u in Eq. ( A27| ) becomes uo = uf for B-orbits. 
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FIGURES 



FIG. 1. Probability density for the natural frequencies {ui u ,lj v ,uj v } at different deformations r] 
in the prolate sodium cluster described in text. 

FIG. 2. M = 0, 1 components of the photoabsorption cross section per valence electron for 
a prolate cluster at different deformations and both in the "single-particle" (shaded black) and 
"RPA-type" (shaded gray) approximations. 

FIG. 3. Fraction of particles on B (full curve) and W (dotted curve) orbits as a function of 
deformation r\ for an oblate sodium cluster. 

FIG. 4. Probability density for the natural frequencies {lo^ , , io^} characteristic of W orbits 
at different deformations r\ in the oblate sodium cluster described in text. 

FIG. 5. As in Fig. |] for natural frequencies characteristic of B orbits. 

FIG. 6. As in Fig. || for an oblate cluster. 

FIG. 7. Photoabsorption cross section per valence electron for prolate and oblate sodium clus- 
ters at several deformations. The left-hand side of the figure shows the results obtained in a 
single-particle approximation (o"(°)). Arrows indicate position of plasmon peak in the spherical 
case. 

FIG. 8. Photoabsorption cross section per valence electron for prolate and oblate sodium 
clusters at several deformations and for N = 25 (shaded black) and N = 254 (shaded gray). 
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